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Abstract 

In this paper we develop a new Peridynamic approach that natnrally includes varying horizon 
sizes and completely solves the ’’ghost force” issue. Therefore, the concept of dual-horizon 
is introduced to consider the unbalanced interactions between the particles with different 
horizon sizes. The present formulation is proved to fulhll both the balances of linear momen¬ 
tum and angular momentum. Neither the ’’partial stress tensor” nor the ’’‘slice” technique 
are needed to ameliorate the ghost force issue in |T]. The consistency of reaction forces is 
naturally fulhlled by a unihed simple formulation. The method can be easily implemented 
to any existing peridynamics code with minimal changes. A simple adaptive rehnement 
procedure is proposed minimizing the computational cost. The method is applied here to 
the three Peridynamic formulations, namely bond based, ordinary state based and non¬ 
ordinary state based Peridynamics. Both two- and three- dimensional examples including 
the Kalthof-Winkler experiment and plate with branching cracks are tested to demonstrate 
the capability of the method in solving wave propagation, fracture and adaptive analysis . 
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1. Introduction 


Peridynamics has recently attracted wide interests for researchers in compntational solid 
mechanics since it provides the possibility to model dynamic fractnre with ease. The crack is 
part of the solution from PD simulation instead of part of the problem and no representation 
of the crack topology is needed. The original PD method was proposed by Silling [2] in 2000 
and has been exploited onwards for extensive applications of mechanical problems including 
impact loading, fragmentation M, composites delamination |S], beam and plate structures 

In PD, the classical balance equations are formulated in an integral form instead of 
partial differential form. The particles interact with each other when the distances between 
the particles are within a threshold value, called horizon. The equations of motion at any 
time t is expressed as 

pu(x, t)= f (u(x', t) — u(x, t), x' — x) — f (u(x, t) — u(x', t), X — x') dDx' + b(x, t) , 

( 1 ) 

where ifx denotes the horizon (spherical domain) belonging to x, u is the displacement 
vector, b denotes the body force, p is mass density in the reference configuration, and f 
is a pairwise force function that computes the force vector (per unit volume squared) [H]. 
Eq. Q shows the key idea of PD to unify continuous and discontinuous media within a single 
consistent set of equations. The governing equation is written in an integral form instead 
of an partial differential form. The crack surfaces are formed as the outcome of motion 
and constitutive models, and there is no entanglement of additional crack kinematics or 
geometry treatment. The fracture behaviour including crack branching and coalescence of 
multiple cracks is captured through the breakage of the bonds between particles. Therefore, 
the need of smoothing of crack surfaces or branching criterion in the extended hnite element 
method (XFEM) [9] , meshless methods [10] or other partition of unity methods (PUM) [TT] 
is completely removed. The extension of PD from 2D to 3D problems is greatly facilitated 
for the computer implementation, which is not always the case in other methods psiini. 


2 


The original PD started with the bond-based formnlation (BB-PD) where the bonds 
behave like springs and independent of each other. BB-PD can be regarded as a special 
case of a more general theory, the state-based peridynamics (SB-PD) [21 [H] which can be 
snited for, theoretically, any type of constitntive model and large deformation analysis. The 
key difference between SB-PD and BB-PD is that in the former, the bond deformation 
depends on collective deformation of other bonds, whereas the bonds in the latter deforms 
independently. The SB-PD was later extended into two types, namely the the ordinary state 
based peridynamics (OSB-PD) and the non-ordinary state based peridynamics (NOSB-PD). 
In all the above types of PD, horizons sizes are commonly reqnired to be constant to avoid 
spnrious wave reflections and ghost forces between particles. However, in many applications, 
the spatial distribntion of the particles with changing horizon sizes is necessary, e.g. adaptive 
refinement, mnltiscale modelling and mnltibody analysis. In other words, in order to achieve 
acceptable accnracy, the entire nnmerical model has to be discretised with respect to the 
highest particle resolntion locally reqnired, and the smallest horizon size used accordingly. 
This is computationally expensive and undesirable. The restriction of horizons being position 
independent practically reduces the efficiency of PD. 

In this paper, we aim to remove the issue of varying horizons and ghost force by de¬ 
veloping a new PD formulation. The new approach is based on the concept of horizon 
and dual-horizon. Though peridynamics has been developed for different types of physical 
fields, e.g. thermal field and fluid field, we confine the present work to solving solid me¬ 
chanics problems. The content of the paper is outlined as follows. §2 begins with stating 
the phenomenon of the ghost forces. In §3, the horizon and dual-horizon are introduced. 
New motion equations with varying horizons are derived based on the dual-horizon concept. 
The balance of linear momentum and the balance of angular momentum of the present PD 
formulation are proved. The applications of the new formulation for BB-PD, OSB-PD and 
NOSB-PD are described in details in §4. In §5, three numerical examples are presented to 
validate the present method. 
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2. Ghost force and spurious wave reflection in peridynamics 


In the original PD theory, the forces exerted on a particle is the snmmation of all the 
pairwise forces from the particles falling inside the horizon of that particle. When the horizon 
sizes are set constant for all particles, the bond force is always pairwise. On the other hand, 
it shonld be noted that the size of the particles that represent the mass qnantity can vary 
[21 [15]. Spurions wave reflections emerge when the horizon sizes vary. We illustrate the 
mechanism of the phenomenon as follows. Consider a particle x' falling inside the horizon of 
particle x, see Figjl} Let fxx' denote the force vector acting on particle x due to particle x', 
where the first subscript x indicates x being the object of force and the second subscript x' 
indicates fxx' being the source of force from x'. Take the Fig. [T]for an example, as x' G iLxj 
the force fx'x ^ 0, x ^ iLx' due to unequal size of horizons, hence fxx' = 0. When computing 
the reactive forces for particle x, the bond forces Fxx' = fxx' — fx'x = —fx'x, which is added 
to Fx. Likewise, when computing forces for particle x', particle x exerts no force on x' as x is 
not inside the horizon of x'. Consequently, the bond force Fxx' only exists unilaterally, which 
is known as the “ghost force” [1] resulting in an unbalanced internal force. The balance of 
linear momentum and balance of angular momentum in this case are violated, and hence 
yields spurious wave reflections in PD simulations. 

Efforts have been made by researchers to explore the possibility of making PD suitable 
for nonuniform spatial discretisation. Bobaru et. al studied the convergence and adaptive 
refinement in ID peridynamics [T6| and multi-scale modeling in 2D BB-PD HZI. Dipasquale 
et al. recently [IB] introduced a trigger based on the damage state of the material for 2D 
refinement of BB-PD. Refined PD was developed in [l9l[T6]|T7l|THl[2n], however it is restricted 
to BB-PD formulation and to the authors knowledge, has only been used for one- and two- 
dimensional problems. In these works, the spurious wave reflection or ghost force problem is 
not solved, and the refinement is performed by checking the spurious reflection is within an 
acceptable range compared to the magnitude of the whole wave. Recently, the partial stress 
was proposed in [1] to remove the ghost force for varying horizons. However, the method 
imposes certain restrictions to the deformation of the body and requires the computation of 
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Figure 1: Force vector in peridynamics with varying horizons. Reactive force —fx'x on x due to x exerted 
on x' as x' S H^. 

partial stress which is complicated. It to certain extent impairs the simplicity of the original 
PD, especially for BB-PD and OSB-PD. Besides, the method does not completely remove 
the ghost force but with a small residual which is believed to be acceptable. Though the 
slice method was devised in the same work, it is applicable only to piecewise constant sizes of 
horizons, and requires additional computation to enforce the consistency between particles 
close to the interface. 

3. Governing equations based on horizon and dual-horizon 

In this section, the concept of horizon and dual-horizon will be introduced to formulate 
the balance equations for particles with varying horizons. It is applied to all peridynamic 
formulations. 

3.1. Horizon and dual-horizon 

We will begin by restating the original concept of horizon in peridynamics. In peridy¬ 
namics theory, the particles interact with each other within a hnite distance. A particle is 
considered to have an influence over other particles within a small domain centering that 
particle. The radius of the domain is known as “horizon”, see Figj^ Particle x (thick solid 

line) is included in the horizons of particles xi, X 2 , X 3 and X 4 (thin solid line), however is 
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not included in the horizons of particles X 5 and xg (dashed line). The concept of horizon 
can be compared to the concept of “nodal support” in meshless methods. 



Figure 2: The schematic diagram for horizon and dual-horizon, all circles above are horizons. The green 
points {xi,X 2 ,X 3 ,X 4 } G whose horizons denote by thin solid line; the red points {xsjXg} ^ ,whose 

horizons denote by dashed line 

Horizon 

The horizon Hy^ is the domain where any particle falling inside will receive the forces exerted 
by X. Hence, x will undertake all the reactive or passive forces from particles in Hy. Hence, 
the horizon can be viewed as passive force horizon. The reactive force acting on x follows 
Newton’s third law. As shown in Figj^ the reactive force —fx'x exerted on x by other 
particles is in opposite direction of the forces applied by x to other particles. 

Dual-horizon 

Dual-horizon is dehned as a union of points whose horizons include x, denoted by H'^ = 
{x' ; X G hfx'}- In the notation of dual-horizon H'y., the superscript prime indicates “dual”, 
and the subscript x denotes the object particle. It can be understood as the set of the 
horizons that belong to the particles who can “see” x in their horizons. As shown in Figj^ 
the dual-horizon with respect to x is the union of xi, X 2 , X 3 and X 4 , whose horizons are 
denoted by thin solid circles. Particles X 5 and xg are not included in the dual-horizon of 
X since their horizons do not include x. In this case, x becomes the object “observed” by 
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the other particles. If x is within the horizon of x', then x' has an active or direct effect on 
X, corresponding to the passive effect in horizon defined previously. Particle x receives the 
active forces from other particles in the dual-horizon of x, and in this sense it is considered 
as “dual” corresponding to horizon. For any point x, the shape of H'^ is arbitrary, and 
depends on the sizes and shapes of horizons as well as the locations of the particles. Note 
that the horizon can take other shapes other than circles or spheres. 

The bond inside the dual-horizon {H'^) is termed as “dual-bond”, and this is correspond¬ 
ing to the bond in the horizon {H^). It can be seen that the bond of one particle can become 
the dual-bond for another particle interacting with it. Note that for each particle, the bond 
and the dual-bond are independent from each other; the same applies to the horizon and 
dual-horizon. It means the bond and dual-bond can break independently in the fracture 
models. For discretisations with varying horizons, often one particle is within the horizons 
of other particles but not vice versa. In this case, a single horizon is not sufficient to define 
the interactions between particles. The concept of two horizons proposed here can solve 
this issue with a simple and direct physical meaning. It naturally takes into account the 
interactions between particles of varying horizon sizes. For models with constant horizons, 
horizon and dual-horizon will degenerate to the horizon in original peridynamics. 

Reaction force by horizon and dual-horizon 

In our dual-horizon peridynamic formulation, the horizons are differentiated between how a 
particle receives and exerts forces with other particles. Under this new concept, computing 
the force fxx' between a pair of particles, denoted as particles x and x', is determined by 
whether x is within the horizon of x', and whether x' is inside the dual-horizon of x, and vice 
versa for fx'x- It can be easily seen that the force density vector has the following property, 

if X G ifx' or x' G H^, fxx' 7^ 0 , 

else fxx' = 0 . (2) 

For any bond between two particles x and x' belonging to a domain denoted as fl in Figj^ 
the direct force fxx' acting on x due to x' can be computed by two approaches as follows. 
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Approach 1 computes the force in terms of x, 


fxx' = 


^0 if x' G i/; 
0 if x' ^ H!^ 
and Approach 2 is formulated with respect to x', 

7 ^ 0 if X G hfx' 


fxx' = 


0 


ifx ^ hfx 



Figure 3: Double summation of force 

For any domain under consideration e.g. the shaded area in Figj^ the computation of 
the direct forces that take place between particles (no reactive force considered yet) shall 
undertake all forces from any particle that belongs to f2. There are two approaches to 
achieve that. The hrst approach is by summing all the forces the particles undertake in the 
dual-horizon, 


E 


^ fxx' Af^x 


,x'em 


ak. 


( 3 ) 


And the second is to add up all the forces of each particle that it applies to other particles, 


E 

x'en 


, xei?^ 


fxx' AFx 


AI4. . 


( 4 ) 
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Since the total force for any fl is independent of the approach chosen to compute it, Eq. ([^ 
and Q shall be equal 

EE fxx' A14'A14 — E E f- AKA 14 ' . ( 5 ) 

xgf2 x'S-ffx x'gfJxSi?,^/ 

Eq. ([^ indicates that changing the summation order from x — )■ x' shall be done together with 
changing from ifx' —t H!^ so as to keep all the variables consistent. When the discretisation 
is sufficiently fine, the summation shall approximate to the integral and Eq. (|^ becomes 


fxx' dl4'dl4 = 


fxx' dl4dl4' 


( 6 ) 


JJJJ 
3 . 2 . Motion equation for peridynamics with horizon-variable 

In the following derivation x' — x will be denoted as and u(x', t) — u(x, t) as 77 , hence 
T] represents the current relative position vector between the particles. The internal 
forces that are exerted at each particle from the other particles should include two parts, 
namely the forces from the horizon and the forces from the dual-horizon. The other forces 
applied to a particle include the body force and the inertia force. Let AI 4 denote the volume 
associates to x. The body force for particle x can be expressed as b(x, f)AI 4 , where b(x, t) 
is the body force density. The inertia is denoted by pu(x, t)AI4, where p is the density 
associated to x. At any time t, for x' in dual-horizon of x, the force vector of fxx' is defined 
as 


fxx' := fxx'(?7, $,) ■ AI 4 ■ AI4' , 


(7) 


where fxx'(?7,^) is the force density in the traditional peridynamics with unit of force per 
volume squared; fxx' is the force vector acting on particle x due to the attraction or repulsion 
from x'. The forces from the dual-horizon are active forces as they are applied to x. The 
total force applied to x from its dual-horizon, H!^, can be computed by. 


fxx' = fxx'(r/,^) ■ AI 4 ■ AI4' . 


( 8 ) 


x'effi 


x'eH' 


Eor any x' inside the horizon of x, the force fx'x acting on x' due to x is defined as 

fx'x := fx'x(-?7, -$.) ■ AI4' ■ AI 4 (9) 
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where fx'x (—'^5 —^) is the force density per volume squared; fx'x is the force vector acting 
on particle x' due to the attraction or repulsion from x ( x' inside horizon of x). The total 
force X exerted on other particles from horizon is the summation of fx'x 


fx'x = fx'x(-r/, -^) ■ AI 4 ■ A14' (10) 

x'er^x x'eHx 

Based on Newton’s third law, the total force x undertaken in i^x takes opposite direction 


- fx'x = - fx'x(-r;, -^) ■ Al^x ■ Al^x' (11) 

x'erfx x'eHx 


By summing over all forces on particle x, including inertia force, body force, active force in 
Eq. (|^ and passive force in Eq. 01 . we obtain the equations of motion 

pu(x,f)A14 = fxx' + [ - ) + b(x, f)A14- (12) 

x'er^i \ x'eHx / 


Substituting Eqs. ([^ and ( pT| into Eq. (12) leads to 


pu(x,f)A14 


^ fxx'(r/,^)AEx'AEx- fx'x(-r;,-^)AEx'AEx + b(x,f)AEx . (13) 

x'eHi x'eHx 


As the volume A14; associated to particle x is independent of the summation, we can elim¬ 


inate AI 4 in Eq. (13), yielding the governing equation based on x: 


pu(x,f)= ^ fx'x(-? 7 , -^)A14'b(x,t) . (14) 

x'er^4 x'eHx 

When the discretisation is sufficiently hne, the summation is approximating to the integra¬ 
tion of the force on the dual-horizon and horizon. 


and 


lim 

AV,,- 5 - 


Y] Uyc'{'n,^)AV^: = [ fxx'(^>^)dlA 
'x'eir; 


(15) 


" x'er/x 


-ri, -^) dl4. 


'er/x 


(16) 
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Thus the integration form of the equation of motion in peridynamics with dual horizon is 
given as 


pu(x,f) = 




fx'x(-? 7 , -^) df 4 ' + b(x, t) 


(17) 


/x'eHi 


'x'eHx 


are 


Eq. (17) is similar to Eq. ([^ of the original peridynamics theory. When the horizons 
set constant, i.e. both horizon and the dual-horizon are equal, the integrations in Eq. 0 
degenerate to the original peridynamics theory. It means the traditional peridynamics can 
be viewed as a special case of the present dual horizon peridynamics. 

About the implementation of the present peridynamic formulation, for any particle x, 
the force density fxx'(?7, in can be determined when calculating the force in for 
particle x'. Therefore, it is not necessary to know exactly the dual-horizon geometry and 
the formulation can be implemented with minor modihcation of any peridynamic codes. 


3.3. Proof of basic physical principles 
3.3.1. Balance of linear momentum 

The internal forces shall satisfy the balance of linear momentum for any bounded body 
f2 given by 


In 


(pu(x,f) - b(x,f))dI4 



n Jx'eH' 



fx'x(-r;, -$,) di4.dK 


n Jx'eH^ 


= 0 . 


(18) 


For simplicities, let fxx' represent fxxK^^^); fx'x the fx'x(——^)- 


Proof: For convenience, we start with the discrete form of the Eq. (14) as 
5^(pu(x,f)-b(x,f))AK = 5^ 5^ fxx.AK'AK-5^ fx.xAK'AK (19) 

xen xeOx'ei^^ xeOx'er/x 

The hrst term on RHS of Eq. (19) can be substituted with Eq. (|^ by changing the 
summation domain and therefore we will get. 


5^(pu(x,f)-b(x,f))AK= fxx'AK'AK-5^ fx.xAKAK , (20) 

xgf2 x'&nxGH^/ xgf2 x'sHx 
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The above transformation is based on the idea that fxx' acting on x is computed by using 
dual-horizon {H'^) of x. And fxx' is calculated by using horizon (-ffx') of x'. Note that in 
both calculations, fxx' remains unchanged and only the dehnition domain changes. Since 
the dummy variables can be relabeled for x <H- x' in the hrst term on RHS, thus the hrst 
term has the same expression as the second term, i.e. 

5^(pu(x,t)-b(x,t))Af4 = 5^ fx'xARxARx'-J^ fx'xARx'Af4 = 0 

As the forces are differentiated here between active and passive force, i.e. an active 
force from x corresponds to a passive force of x', the forces is always pairwise but with 
opposite direction and of the same magnitude. Hence, it is natural that the balance of 
linear momentum is satished. 

3.3.2. Balance of angular momentum 

Let y be the deformation vector state held dehned by 

y[x,t](^) = y(x,f)-y(x,t) Vx G Li, ^ = x - x, x G iLxA > 0 (21) 


u(x,t) = y(x,t) -X , 


( 22 ) 


where y(x, t) refers to the current conhguration coordinate for x in material conhguration, 
u(x, f) is the displacement for x, is the horizon. Thus y (x' — x) is the image of the bond 
x' — X under the deformation. 

To satisfy the balance of angular momentum for any bounded body 17,it is required that: 


y X (pu(x, t) - b(x, t))dHx 

= / yx f / fxx'(r/,^)dl4'- [ fx'x(-? 7 ,-^) dl4') df4 

Jn / 

y X fxx'(^7, i) df4'dl4 - [ [ y X fx'x(-^7, -$,) dlA'dlA 



n Jx'&H’ 


n Jx'eH^ 


= 0 


For simplicities, let fxx' represent fxx'(^,^), and fx'x present fx'x(——^)- 
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( 23 ) 



Proposition: In the dual-horizon peridynamics, suppose a constitutive model of the form 


f = f(y,A) (24) 

where f ; V —)■ V is bounded and Riemann-integrable on H and V is the vector state; A 
denotes all variables other than the current deformation vector state that may depend on 
for some particular material. 

If 


or in discrete form 



X fx/x dI4' 


0 


VyeV, 


(25) 


y(x' - x) X fx'x AI4' = 0 Vy G V , (26) 

x'erfx 


where fx'x = fx'x(u(x', f) — u(x, f), x'—x) , the force vector density acting on x', and is the 


horizon of x. Then, the balance of angular momentum, Eq. (23) holds for any deformation 
of for any given constitutive model. 

Proof: As the summation over the domain 17 is equivalent to the integrand of the same 
expression over that domain, we use the discrete form as 


y X (pu(x, t) - b(x, f))AI4 
= y^(x + u) X (pu(x, t) — b(x, t))AVx 

x60 

= EE(- + u) X fxx' AI4'AI4 - EE (x -h u) X fx'x AI4'AI4 (27) 

xeOx'eHx xsOx'eHx 


The hrst term on the RHS of Eq. (27) can be rewritten according to Eq. (|^ by changing 
the summation domain and thus it becomes 


y X (pu(x, t) - b(x, f)) AI4 

xgf2 

= E E(- -I- u) X fxx' AI4;/AI4 — EE (x -h u) X fx'x AI4'AI4 • (28) 

x'gf2 xSf2 x'SiTx 
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Swapping the dummy variables for the hrst term on RHS between x t-)- x' in the double 
summation, then the hrst term takes the same summation to the second term on RHS 

^ y X (pu(x, t) - b(x, t)) AVk 

xgf2 

= EE( x' + u') -EE(- + u) X fx'x A14;/AI 4 

= ^ ((x' + u') - (x + u)) X fx/x ARxO^Kc 

xGil x^Gffx 

= (y' - y) X fx^x ARxOARx 

xefJ x'ei?x 


= ^ y (x' - x) X fx'x A14') ARx = 0 

xefJ x'ei?x 

or in an integration form for sufficiently hne discretisation as 


(29) 


y X (pu(x, t) - b(x, i))dl4 ^ 0 


(30) 


' xef2 


Therefore, the angular momentum over the entire analysis domain is satished. □ 



BB-PD andOSB-PD 



NOSB-PD 


Figure 4: force vector of BB-PD,OSB-PD and NOSB-PD 


It can be seen that the conservation of angular somehow depends only on the horizon; 
the dual-horizon is not involved. The latter is only needed in the Eq. (14) and Eq. (0. 


The bond-based, ordinary based and non-ordinary based peridynamics all satisfy the angular 
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momentum in the original horizon concept (see proof in [IS]). This conclusion can be also 
illustrated for the BB-PD and OS-PD since the internal forces fxx' and fx'x are parallel to 
the bond vector in the current configuration, see Figj^ 


4. Dual-horizon peridynamics 


The dual-horizon formulation will now be applied to all existing peridynamics, namely 
BB-PD, OSB-PD and NOSB-PD. The calibration of constitutive parameters with respect 
to the continuum model and some issues concerning the implementations will be discussed 
for 2D and 3D problems. 


4-.1. Dual-horizon bond based peridynamics 

In the bond-based peridynamics theory mi, the pair force function is expressed as 

= (31) 

For a micro-elastic homogenous and isotropic material, f(77,^) can be specified as 

f(T/,^) = C5- 11 ^ I ^11 1 (32) 

where c is the micro-modulus, and s is the bond stretch calculated by 



The energy per unit volume in the body at a given time t is given by [S] 


^ = 11 (34) 

J Hx 

For the BB-PD theory, by enforcing the strain energy density being equal to the strain 
energy density in the classical theory of elasticity la El, and by setting the influence function 
(udl^ll) = 1, we obtain 


(^(S\ _ __ 

^ ^ “ tt 6%1 + iy){l - 2iy) 

r'l'A'l — 3i^ 

“ 7r54(i - 2z/) 


plane stress 
plane strain 

3D . (35) 
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Note that the value of C'((5) takes half of the micro-modulus c used in the horizon-constant 
BB-PD since the bond energy for varying horizon is determined by both horizon and dual¬ 
horizon. It can be easily seen that when the horizon takes the same value as the dual-horizon, 
the bond energy between two particles is reduced to that of the original BB-PD. 

Let Wo(^) = C{5)sl{5)^/2 denote the work required to break a single bond, where So(h) 
is the critical bond stretch. By breaking half of all the bonds connected to a given particle 
along the fracture surface and equalizing the breaking bonds energy with the energy release 
rate Gq [mil], we can get the expression between the energy release rate Gq and the critical 
bond stretch So(5): 


SoW = 1 

UttGo 
\ l 9E6 

plane stress 

SoW = 1 

1 birGo 
\l 12E6 

plane strain 

SoW = 1 

/5Go 
\l 6E6 

3D . 


(36) 


Both the micro-modulus ^((f) and the critical stretch so(^) are derived from the local con¬ 
tinuum mechanics theory, and they depend on the horizon radii for variable horizons. 

In the implementation of the bond-based peridynamics, fracture is introduced by remov¬ 
ing particles from the neighbour list once the bond stretch exceeds the critical bond stretch 
So- In order to specify whether a bond is broken or not, a history-dependent scalar valued 
function /i is introduced [8|, 





if s(t', < So for all 0 < < f, 

otherwise. 


(37) 


The local damage at x is dehned as 

i).(x,i) = 1 - 


Ih. dU 


(38) 


The damage formulation of Eq. (38) is also applicable to OSB-PD. For any particle with 
dual-horizon and horizon (Lfx), the active force fxx' and the passive force fx'x in Eq. 
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(14) or ( [ItI ) are computed by the following expressions, respectively 

r/ + ^ 


fxx' - C'(^x') ■ "Sxx' 

fx'x = C'(^x) ■ Sxx' ■ 


N+^ir 

-(y/ + ^) 

l^ + ^ll ’ 


Vx' G Hi 


Vx' e , 


(39) 

(40) 


where (^(hx') and C{6x) are the micro-modulus based on hx' and ^x computed from Eq. (35) 
respectively, and Sxx' is the stretch between particles x and x'. 

4-2. Dual-horizon ordinary state based peridynamics 

The concept of “state” for peridynamics was firstly introduced in [15]. A state of order 
m is a function 

(41) 

where H is the horizon domain, denotes the set of all tensors of order m, (□) indicate 
the vector to which a state operates. In this paper, the scalar state and vector state are use 
if not otherwise specified. 

The 2D OSB-PD formulation of bond force can be derived by following the similar 
procedure as that was described in 3D. Both 2D and 3D OSB-PD can be written in a 
unihed expression taking into account the dimension number as 


nKO ... ^ nin + 2)G 


(42) 


m m 

where n G {2, 3} is the dimensional number, K is the bulk modulus, G the shear modulus, 

i = e(^) - —, m = 6 = — f n;(^) ^ • e(^)dV'^, e(^) = 

^ ^ Jh 

11^ + ^11 ~ 11^II dV^ are the area in 2D, and volume in 3D, respectively. For any particle 


X, fxx' and fx'x in Eqs. (14) and (17) can be calculated by 

r} + ^ 


and 


fxx' = t{^) 


fx'x = t'{-0 ■ 


i^+^ir 


Vx' G Hi , 


Vx' G ifx 


(43) 


(44) 


l^ + ^ll 

Here ^ = x' — x, and x and x' are the coordinate vectors for x and particle x' in the material 
conhguration, respectively. 
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4-3. Dual-horizon non-ordinary state based peridynamics 

NOSB-PD uses the deformation gradient from classic continuum mechanics. It offers the 
possibility to consistently include existing constitutive models. With the use of the shape 
tensor proposed in the NOSB-PD, the fundamental tensors in continuum mechanics can 
be easily introduced in peridynamics, including the deformation gradient tensor or velocity 
gradient |22] . The shape tensor for particle x is dehned as 


Kx=/' (45) 

^x'e/Tx 

where ^ = x' — x is the bond in the reference conhguration, is the horizon for particle 
X and is the influence function. The deformation tensor for particle x is given by 

F, = ^ / w(i)y{i) ® ? dF,, . K;‘ , (46) 

where y := y(x, t) are the spatial coordinates, x are the material coordinates, and ^ = x'—x. 
The spatial velocity gradient Lx for particle x is dehned in the current conhguration using 
the chain rule. 


dy dx dy ^ ’ 


(47) 


dy, 


where v := v(x, t) = —(x, f) is the velocity vector and Fx is the rate of deformation 
gradient. The non-local form of the velocity gradient can be written as 


Lx; = 


w(^)v(^) dFx' 


/x'eHx 


a;(^)y(^) ® ^ dI4/ ■ K 


-1 


'x'eHx 


[ a;(^)v(^) (g)^ dI4/ • Kx f / c^(^) y(^) ® ^ dI4 

V./x'eHx 


-1 


a;(^) v(^) g)^ dFx 


'x'eHx 


^^(^)y(^) di4 


'x'eHx 


-1 


(48) 


where ^ = x' — x, v(^) = v(x',f) — v(x, f) and y(^) = y(x',f) — y(x, f). It can been seen 
that the velocity gradient does not require shape tensor, which means the shape tensor is 
not necessary to dehne velocity gradient. The incremental deformation gradient can also be 
obtained without shape tensor K. 
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Let Au denotes the displacement increment with respect to the previous time step. The 


incremental spatial deformation gradient Cx in continuum mechanics is defined as 

dy dy 5x 

The non-local counterpart is given by 

Cx:= / n;(OAu(^)0^dl4'-(/ (^) ® ^ dV;,)"' • 

Many material models in non-ordinary state based peridynamic accounting for the geomet¬ 
rical and material non-linear problems are solved on the incremental spatial deformation 
gradient. The non-ordinary bond force is computed based on the Piola-Kirchhoff as 

t,,, = u,(OPx'K;‘-« (50) 



Px = Jx = det Fx 


(51) 


where ^ = x' — x,(Jx is Cauchy stress tensor. 

Let us dehne the deformation gradient 

F:, := Fx • Kx = [ u{^)y{0 ® ^ dFxi ■ Kx = [ u{^)y{0 ® ^ . (52) 


Therefore 


Jx = detFx = ^^ 
det Kx 


By substituting Eqs. (53), (52) and (51) to (50), we obtain 


detFx 

detFx 

= <mk; 

det Fx 
det Kx 


= ij(?) 


^■K-A^ 



■^x-(FA 

Kx') 

■ K-^ 

■ cTx ■ f:,-^ 

Kx 

KxA^ 

■ ■ F-^ 

■ 



(53) 


(54) 
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The shape tensor is only needed for = detFx, while the compntation of the bond force 
density does not reqnire the shape tensor. For any particle with horizon ifx and dnal- 


horizon H'^ , the active force fxx' and the passive force fx'x in motion Eq. 
be calcnlated by 


(14) or (17) can 


fxx' - Txx' 


Vx' G H' 


(55) 


and 


fx'x = Dx'x Vx' G FTx , (56) 

where ^ = x' — x, and x and x' are the coordinate vectors for particle x and particle x' in 
the material conhgnration, respectively. 

The snmmary of three kinds of dnal-horizon peridynamics see Tablej^ 


Model 

original Peridynamics 
/y^(fxx' - fx'x)dl4' 

dnal-horizon Peridynamics 
jiji fxx'dV^' Jrtx fx'xdV^' 

BB-PD 

fxx^ fx^x C * 5xx^ ' ^^5 ^ 

^ + ^ 

fxx' = C'(^x') • Sxx' ■ n, Vx' G 

fx'x = C{5^) ■ Sxx' ■ (-n),Vx' G i7x 



V + ^ 


VV lit71 Kj 11 . 

ri + i 

OSB-PD 

fxx' = t{i) ■ n, 

fx'x = ■ (-n), Vx' G i7x 

where n = -— 

r] + i 

fxx' = ^(^) ■ n,Vx' G 

fx'x = ■ (-n), Vx' G TTx 

where n = -— 

r] + i 

NOSB-PD 

fxx' = Txx') 

fx'x = Tx'x) Vx' G i7x 

fxx' = Txx',Vx'G77;; 

fx'x = Tx'x)Vx' G 


Table 1: Comparison of the original Peridynamics and dual-horizon Peridynamics 
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5. Numerical Examples 

5.1. Two-dimensional wave reflection in a rectangular plate 

Consider a rectangular plate with dimensions of 0.1x0.04 (see Figj^. The Young’s 
modulus, density and the Poisson’s ratio used for the plate are E = 1, p = 1 and 1 ^ = 0, 
respectively.Note that this is 1-D model. The initial displacement applied to the plate is 
described by the following equation 


X 


uo{x,y) = 0.0002exp[-(^^)^], vo{x,y) = 0, x e [0,0.1],?/ e [-0.02,0.02] , (57) 


where uq and vq denote the displacement in the x and y directions respectively. The wave 
speed is n = \/E/p = 1 m/s. At any time step, the L 2 error in the displacement is given by 


u — u 


|err||^^ = 


analytic | 


U 


analytic | 


(58) 


with 


u = 


U ■ U dflr 


Oq 


Three models for solving this problem, namely model A, B and C, are devised,see Figs,6(a] 
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Figure 5: setup of the plate 
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and 6(b) In model A, the plate is discretised with 4000 particles and all particles have the 


same horizon sizes of 0.001 m, see Fig,6(a) For models B and C, we adopted a discretisation 


of varying horizons as shown in Fig, 6(b) Verlet integration method is adopted.The horizon 
radins associated to each particle is set as 3 times the particle sizes for all models. Therefore, 
along the interface between the coarse and fine discretization, the horizon sizes vary. The 
minimal particle size of 3 models is Ax =5e-4m, yielding a critical time increment Atmax = 
Ax/n=5e-4s. The physical compntational time step is set 0.5 second with the time increment 


AA„,„=5e-4s. 


As BB-PD can not simnlate the model with Poisson ratio z/ = 0, the OSB-PD model is 
adopted. Model A was solved with the OSB-PD with constant horizons. Model B is the 
original OSB-PD with variable horizon (withont additional treatment for ghost force), and 
model C is onr dnal-horizon formnlation for OSB-PD with variable horizon. The parameters 
nsed for the three models are listed in TablelH 


Model 

Ax = Ay 

0 

CO 

Particle nnmbers 

A 

0.001 

3 

4000 

B 

0.001,0.0005 

3/1.5 

8628 

C 

0.001,0.0005 

3/1.5 

8628 


Table 2: Values of the Peridynamic parameters for the mixed model 


Fig j7|[8|^ show the displacement in the x-direction along the x-coordinate of all particles. 
The red star points are particles with Ax = 0.0005m, the bine dot particles refer to Ax = 
0.001m. The displacement in model C (present dnal-horizon PD) is almost identical to 
that of model A (PD with constant horizon),as shown in Figj^ and Figj^ Spnrions wave 
reflections are observed for Model B, as shown in Figj^ At step 650, several wave peaks are 
observed in Model B and the maximnm displacement is lower than that of Models A and 
C. Therefore, the displacement wave in Model B were affected by spnrions wave reflection 
and the resnlts deviate greatly from resnlts of A and C. 

We compnte the L 2 errors of the wave prohle along the line-AD (see Fig. |^. The L 2 
error (see Table. of the present peridynamic formnlation, can achieve the almost the same 
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(a) Particles distribution of model A 



(b) Particles distribution of model B and C 


Figure 6: The discretizations of model A,B and C. 
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X displacement (10 m) 



X coordinate (m) 


Figure 7: Displacement wave versus x coordinate of model A at step 650 
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X displacement (10 m) 



0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 

X coordinate (m) 


Figure 8: Displacement wave Ux versus x coordinate of model B at step 650 
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X displacement (10 m) 


1.2 


• A x=0.001 m 
+ Ax=0.0005 m 
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X coordinate (m) 


Figure 9: Displacement wave Ux versus x coordinate of model C at step 650 
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accuracy as the conventional peridynamics with constant horizon. 


Model 

step 100 

step 200 

step 650 

step 990 

A 

0.0410 

0.0452 

0.1277 

0.2498 

B 

0.1149 

0.1738 

0.7638 

1.3233 

C 

0.0447 

0.0500 

0.1253 

0.2871 


Table 3: L 2 error of displacement for the initial condition (Gauss distribution of displacement) before and 
after the wave reflection 

Displacement and velocity of monitor points 

Two monitor points (A,B) are selected to evaluate the spurious wave reflection of Models 
B and C as shown in Fig. The displacement and velocity curves show the spurious wave 
reflection exists in model B while not exists in model C, see FigJTOj 

5.2. Kalthof-Winkler experiment 

To study the performance of the present formulation for fracture modeling, the Kalthof- 
Winkler experiment is exploited. The geometry and model used for the test is depicted in 
FigJTH the thickness of the specimen is set as 0.01m. In the experiment, the evolution of 
crack pattern was observed to be dependent on the impact loading velocity. For a plate made 
of steel 18Nil900 subjected to an impact loading at the speed of Vo=32 m/s, brittle fracture 
was observed [ 23 ]. The crack propagates from the end of the initial crack at an angle around 
70° vs. the initial crack direction. In [23], the BB-PD with constant horizon was used to 
test this example. For convenience of comparison, the present dual-horizon formulation is 
also applied to the bond-based peridynamics to test the example. The material parameters 
used are the same as in [23], i.e. the elastic modulus E = 190GPa,p = 7800kg/m^,z/ = 0.25 
and the energy release rate Go = 6.9e4 J/m^. The impact loading was imposed by applying 
an initial velocity at Uq = 22 m/s to the first three layers of particles in the domain as shown 
in FigJTT]all other boundaries are free. The plate is discretized with two different particle 
sizes, namely Axcoarse = 1.5625e-3m for the coarse subdomain and Axdense = 0.5Axcoarse = 
7.8125e-4m for the fine subdomain located in the left down corner of the model, see FigJTTj 
Along the thickness of the plate (in direction perpendicular to the plane surface), four layers 
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(c) X velocity of model C 


(d) y velocity of model C 



Time(s) 



Time(s) 


(e) X displacement of model B 


(f) y displacement of model B 



Time(s) 



Time(s) 


(g) X velocity of model B 


(h) y velocity of model B 


Figure 10: The 


displacement and velocity curve of monitor points A and B for model B and C 
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of particles in the coarse subdomain and eight layers in the hne subdomain are employed. 
The total number of particles is 57968. The crack propagation speed is computed by 
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Figure 11: Kalthof-Winkler’s experimental setup 


Vi-o.b — 


Xz - Xz_i 


(59) 


tl — tz-l 

where x; and xz_i are the positions of the crack tip at the times ti and tz_i respectively. An 
expression for the Rayleigh wave speed cr [2l] is given as 


c. 


0.87 + 1.12z/ 


(60) 


where u is the Poisson’s ration, Cg = s/JjJp is the shear wave speed and /i is the shear 
modulus. The crack starts to propagate at 26.3 /is. The highest crack speed reached is 
1530 m/s, about 54.4% of the Rayleigh speed(2799.2 m/s). The average crack speed is 
1077 m/s for the hne subdomain and 1094 m/s for the coarse subdomain. During the crack 
propagation, the crack paths in the hne and coarse subdomains are nearly symmetrical, as 
shown in Figs. 1^, 1^ and M It can be seen from FigjT^that the crack speed is very close 
to that predicted by the PD formulation using constant horizons (Aa^uniform = 1.5625e-3m). 
The crack propagation initiated at an angle of 65.7° in the hne subdomain with respect to 

the original crack and 65.8° in hne subdomain. Therefore, it can be concluded that with the 
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present formulation, the transition from fine to coarse of horizons has almost no influence 
on the crack propagation speed. 



Figure 12: The crack pattern of Kalthof-Winkler plate by the present dual horizon PD at step 350 



Figure 13: The crack pattern of Kalthof-Winkler simulation by the present dual horizon PD at step 650 


5. 3. Adaptively refined peridynamics 

The present dual-horizon Peridynamic formulation provides the possibility for adaptive- 
rehnement within a simple and unihed framework. Two examples of adaptively rehned 


peridynamics will be tested in this section: the Kalthof Winkler experiment in section 5.2 


and plate with pre-crack subjected to traction. The threshold for the adaptive rehnement 
is determined by both the damage-state criteria and energy state, as proposed in [TH]. The 
adaptive rehnement procedure consists of two steps. 

(1) Search for the particles that exceed the threshhold values. Note that the rehnement is 
not only applied to the particles above the threshold but also to the neighbouring particles. 
In this way, it can be ensure that the crack tip always remains inside the rehned zone. 
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Figure 14: The crack pattern of Kalthof-Winkler simulation by the present dual horizon PD at step 875 



Figure 15: The crack speed of Kalthof-Winkler simulation by the present dual horizon PD 
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(2) Split the particle, named as parent particle, into small particles, named as child particles. 
The properties that will be mapped from the parent particles to the child particles include 
the mass, volume, coordinate, displacement and velocity. 

The method to split particles for structured discretization is shown in Fig. Splitting the 
particle results in halving the maximum time step. In all examples, we restrict each particle 
being allowed to split only once in the entire analysis. 



parent particle child particles 



Figure 16: Particle splitting for Adaptive refined peridynamics 


5.3.1. 3D adaptively refined Kalthof-Winkler simulation 


We now repeat the Kalthof Winkler experiment from section 5^ However, the initial 
discretization is based on uniform horizon size. The initial particle size is Axjnitiai =1.5625e- 
3 m. After refinement, the particle size is reduced to Axrefined =7.8125e-4 m around the 
crack. The total particles is increased from 32,768 to 54,860 at the end of the simulation. A 
uniform rehnement would result 32,768x8=262,144 partilces. 

The crack patterns at certain steps are plotted in Figsj^ and W The crack tip is 
always contained inside the refined zone. According to Fig. the maximum crack propa¬ 
gation speed reaches 1268 m/s and the average speed is 1047.6 m/s. The crack propagation 
takes place at an angle of approximately 66.5° with respect to the initial crack direction. 
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Figure 17: The crack pattern of adaptive refined Kalthof-Winkler experiment at step 420 



Figure 18: The crack pattern of adaptive refined Kalthof-Winkler experiment at step 600 



Figure 19: The crack pattern of adaptive refined Kalthof-Winkler experiment at step 925 
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Figure 20: The crack speed of adaptive refined Kalthof-Winkler experiment 

5.3.2. Plate with pre-crack subjected to traction 

In the last example, we will test the capability of the new formulation in modelling the 
crack branching. Hence, a plate with pre-crack subject to traction is considered as shown 
in Figj^ The traction remains constant during the entire simulation. The same example 
has been well studied in |2S] with BB-PD and also in [2B] using the XFEM. 

The material parameters of the plate (soda-lime glass) are E = 72 GPa, p = 2440 kg/m^ 
and energy release rate being Gq = 135 J/m^ [25]. Two models, namely Case 1 and Case 2, 
are set up for comparison. Case 1 is solved by using the traditional BB-PD with constant 
horizons, while Case 2 is modeled by an adaptively rehned BB-PD using the present dual¬ 
horizon formulation. The plate is assumed under plane stress condition and all simulations 
are carried out in 2D. The particle sizes chosen are 5 x 10“^ m for Case 1 and 1 x 10“^ 
m for Case 2, respectively. Case 2 uses an adaptively rehned model, and once a particle is 
rehned, the minimal particle size will become the same as in Case 1. The particle number 
are 16,000 for Case 1 and 4000-6424 for Case 2, respectively. The Rayleigh speed for the 
material is 3102 m/s according to Eq. ( |60| ). The max crack speed is 1881.4 m/s (60.6% of 
cr) for Case 1 and 2184.1 m/s (70.4% of cr) for Case 2, respectively. One possible reason 
for the diherent max crack speeds is the deviation of the mapping method. 

For Case 1, the crack starts to propagate at 11.9 /rs, and the hrst crack branching point 
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B\ occurs at 24.5 /rs at the speed of 1043.6 m/s and the second crack branch point B2 takes 
place at 40.9 /rs at the speed of 1147.1 m/s. It is observed once branching initiates, the crack 
propagation speed decreases. This can be explained that the energy release to create more 
fracture surfaces decelerates the propagation speed. Afterward, the crack speed increases. 
For Case 2, the initial crack started to propagate at 10.3 /rs, first branching point B1 at 
19.4 /rs with 1247.6 m/s, and second branching point B2 at 34.2 /rs with 1330.6 m/s. The 
crack pattern is similar compared with Case 1. 



Figure 21: Setup of the pre-cracked plate under traction load 
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Figure 22: The crack pattern of pre-cracked plate with uniform horizons 

6. Conclusions 

This paper contributes to the development of peridynamics with varying horizons. There¬ 
fore, the interactions between particles are based on two independent horizons, passive force 
from the horizon and the active force from the dual-horizon. The spurious wave reflection 
and ghost force in the conventional peridynamics is naturally eliminated. Based on the 
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Figure 23: The crack pattern of adaptive refined pre-cracked plate 



Figure 24: The crack speed of pre-cracked plate under traction 
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new concepts of horizons, the motion equation of dual-horizon peridynamics is derived. We 
show that the balance of linear momentum and angular momentum is satished. The key 
difference of the present peridynamics formulation from the conventional one is the way of 
computing reactions forces. Hence, the dual-horizon peridynamics can be implemented in 
any existing peridynamics code with minimal efforts. Three numerical examples show the 
present method is free from spurious wave reflection, and the accuracy is retained along the 
interface where horizon sizes undergo sudden changes. The method also shows its capability 
for fracture problems including crack branching. A simple h-adaptive scheme is proposed 
to improve crack paths resolution, both for 2D and 3D case. We also intend to develop 
efficient error estimate for adaptivity to drive the adaptive rehnement. The present method 
also shows the potential in multiscale analysis where the models at different length scales 
can be bridged by using different horizon settings. 
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